HOMEWORK 2 - PROGRAMMING PORTION¶

Genevieve Mortensen 2023-03-18

Our purpose to programmatically investigate mouse brain cortex scATAC-seq data using methods learned in class.

scATAC-seq stands for single cell Assay of Transposase Accessible Chromatin sequencing. This type of sequencing attempts to capture epigenomic information at a single-cell resolution in an efficient way. The method involves Tn5 transposase nicking open chromatin and ligating adapters to the fragmented chromatin. The resulting library is amplified via PCR and then next-generation sequenced.

scATAC-seq data is known to be sparse and noisy, making data analysis a challenge. So, let's practice!

1.1) Download data¶

Data is downloaded from the following GEO website and unzipped using 7Zip File Manager: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126074

The barcodes correspond to the cells and the peaks correspond to the regions of open chromatin, which can help identify regulatory domains.

Open chromatin indicates that the cell was using the information encoded in that region because it was exposed to the transcriptional and translational cell machinery, indicating some epigenetic importance.

Chromatin peaks with high counts for a cell barcode mean that the region was highly expressed.

            For our dxn matrix, we have d = 244544 features and n = 10309 samples

We will need to transpose the data accordingly and use our nxd matrix when we perform our analysis using EpiScanpy functions because of the way the Anndata class structures data. Please see HW1 for the in-depth explanation.

1.2) Analyze data using EpiScanpy¶

Following the Buenrostro tutorial and the EpiScanpy API (respectively referenced below), we will unpack the data and visualize findings.

https://nbviewer.org/github/colomemaria/epiScanpy/blob/master/docs/tutorials/Buenrostro_PBMC_data_processing.ipynb

https://episcanpy.readthedocs.io/en/latest/api/index.html

We don't have metadata associated with this group of data, so we can move on to preprocessing. Keep in mind that any analysis covered in the tutorial which must use metadata will not be performed in this experimental analysis.

Preprocessing¶

Quality Control¶

The data looks okay in terms of noise and distribution. We need to set a minimum number of cells to keep and features as well.

Let's look at the distributions which are created when we hypothetically remove some minima of features and cells.

For our features, we'd want a normally distributed plot capturing most of our data. For our cells, we'd expect a right-skewed graph capturing most of our data.

Looks like we can filter at our experimentally found minima.

Now, we need to identify the most variable features.

EpiScanpy will give us the a plot showing a variability score for our features. We need to select a cutoff after the elbow of the distribution oof variability scores.

Great, we have chosen cutoff's just around the elbows. Now we can save the raw data and preprocess the data such that we only retain the highly variable features for dimension reduction and cluster identification.

Our data looks to have no outliers from the central distributions, so we don't need to filter anymore. Let's move on!

Normalization¶

Before normalizing, we can look at the binary matrix UMAP. We need to compute the PCA, UMAP neighbors, and TSNE embeddings.

We don't have metadata associated with this dataset, so we cannot visualize cell types like the tutorial does.

            Our matrix now has 10202 cells and 60213 features

Great! Now let's visualize everything

Visualization¶

Let's take a look at our normalized and dimensionally reduced data.

The first plot shows that most of the variation retained by the first two principal components, but clusters poorly.

The second plot shows rankings of each feature to the variation of each principal component.

The third plot shows the rankings of each principal component to the variation of the dataset.

Now that we've checked out the variation in the non-normalized data, let's do dimension reduction on the normalized data and visualize that.

Now we are getting somewhere! We see better clustering of data.

We still don't see great clusters, though, and the number of features recovered per cell vary significantly. Let's log transform to help massage the data into something more workable for identifying clusters, then perform dimension reduction again, then plot the data again.

Clustering¶

Let's do four different kinds of clustering and then compare them.

We need to look at the resolution we need to obtain a certain number of clusters. At a resolution >= 1, we get 7 clusters. We want to retain a reasonable number of clusters. We will use this to set the number of clusters for our kmeans and heirarchal clustering algorithm plotting, which require a positional argument about the number of clusters.

Louvain clustering appears to be the most distinct, but we should run some quantitative metrics to be sure. We don't have metadata associated with the clustering, so it is possible that the following results may be conflicting.

Comparison of clustering¶

We need to check to see which clustering algorithm was the best. First, we will use Adjusted Rand Index, which computes the similiarty of clusters by comparing cells which are assigned to the same cluster. Then, we will use a Homogeneity score, which scores more highly for clustering algorithms where cells have the same ground truth labels. Lastly, we will compute the Adjusted Mutual Information, which measures similarity between two labels of data for the same information while accounting for chance.

For all 3 methods, the closer the score is to 1, the better the quality of the clustering.

            Quantitatively, it appears heirarchal clustering was the best at clustering the features.

            Qualitatively, it appears Louvain clustering was the best. 

            What do we trust?

These conflicting results are because we don't have any metadata associated with our dataset, so the clustering is only dependent on the raw features and not on some relevant categorization. Therefore, the metrics we calculated are actually not as informative as they could be.

This naive approach isn't helpful in a biological sense, so we can use our qualitative skills to perform more analysis and make sense of our data.

Ultimately, it appears we have two clusters- 0 and 6- which have similar features (e.g. chr9: genomic-coordinates) as each other with similar contributions to the characterization of those clusters, meaning they are indistinct clusters. This indicates our data could be clustered better because these two clusters are so similar. Converesely, cluster 1 is not well characterized by any of the features which comprise it, showing that these features only belong in this cluster because they could not belong in any other, indicating a cluster where the relation of features is weak.

On the other hand, clusters 2, 3, 4, and 5 appear distinct and well characterized by their respective features given their high ranking. This means that we can rely on the Louvain algorithm's distinction between the clusters and therefore the features, indicating the regions are all associated with distinctly different regulatory purposes within the cell.